    
    logit <- function(prob){ log(prob) - log(1-prob)}
    
    expit <- function(logodds){ 1/(1+exp(-logodds))}
    
    
    getlogop <- function(p0,p1){
        log(p0) + log(p1) - log(1-p0) - log(1-p1)
    }
    
    getlogrr <- function(p0,p1){
        log(p1) - log(p0)
    }
    
    MyAttach = function(list){  # attach a list to the global environment
        for(i in 1:length(names(list))){
            assign(names(list)[i], list[[i]], envir=.GlobalEnv)
        }
    }
    
    
    
    NameFirst = function(names){
        
        sapply(names,  function(name) strsplit(name,split=".",
                                               fixed=TRUE)[[1]][1])
        
    }
    
    ParaValues = function(paras){
        
        ###     Create an array containing all the parameter value combinations
        para.string = paras[[1]]
        for(i in 2:length(paras))   para.string = outer(para.string,paras[[i]],"paste")
        
        ###     Extract numerical values from strings
        Extract = function(name) as.numeric(strsplit(name,split=" ",fixed=TRUE)[[1]])
        para.value = apply(para.string,1:length(dim(para.string)),Extract)          # The first dimension contains the results
        dimnames(para.value)[[1]] = NameFirst(names(paras))
        dim.order  = c(2:length(dim(para.value)),1)
        return(aperm(para.value, dim.order))
        
    }
    
    
    
    ## Function for checking if two things are equal within numerical precision
    same <- function(x,y,tolerance=.Machine$double.eps^0.5){abs(x-y)<tolerance}
    
   
    
    
    ## Print error message:
    PrintErrorMessage = function(opt){
        if((opt$convergence>0)){
            print(paste("Convergence not achieved",opt$message,
                        "  code:",opt$convergence))
        }    
    }
    
    
    Optimize = function(startpars,objective,Hessian,max.step = 5000){
        # input
        # startpars:          parameters to start with
        # neg.log.likelihood: objective function to minimize
        # Hessian:            logical; whether or not calculate the Hessian
        
        opt <- opt2 <- optim(startpars,objective,hessian=Hessian)
      # opt <- opt2 <- optim(startpars,objective,hessian=Hessian, control=list(trace=1))
        if(opt$convergence != 0){
            cat("The default method does not converge in 500 steps. \n")
            
            diff <- 1
            while(diff > 1e-6){
                startpars <- opt2$par
                opt   <- optim(startpars,objective,control=list(maxit=100000))
                set.seed(opt$counts[1])
                startpars <- opt$par
                opt2  <- optim(startpars,objective, method="SANN",
                               control=list(maxit=max.step))
                
                diff  <- (opt$value - opt2$value)/opt$value
                cat("We compare the default method with the SANN option, the relative difference is ", diff, "after ",opt$counts[1],"steps \n")
            }
            
            
        }  
        
        return(opt2)
    }
    
    
    
    OptimizeRanStart = function(startpars,objective,max.step){
    # For comparison only
    # Not for actual implementation
        
        rep.time = 10
        results  = matrix(,rep.time,length(startpars)+3)
        
        for(i in 1:rep.time){
            set.seed(i)
            startpars <- rnorm(length(startpars),0,1)
            while(!is.finite(objective(startpars))){
                startpars <- rnorm(length(startpars),0,1)
            }
            opt   <-  Optimize(startpars,objective,FALSE,max.step)
            results[i,] = c(opt$par,opt$value,opt$counts[1],opt$convergence)
        }
        
        
        return(results)
    }
    
    OptimizeRanStartDefault = function(startpars,objective){
        # For comparison only
        # Not for actual implementation
        
        rep.time = 10
        results  = matrix(,rep.time,length(startpars)+3)
        
        for(i in 1:rep.time){
            set.seed(i)
            startpars <- rnorm(length(startpars),0,1)
            while(!is.finite(objective(startpars))){
                startpars <- rnorm(length(startpars),0,1)
            }
            opt   <- optim(startpars,objective,control=list(maxit=100000))
            results[i,] = c(opt$par,opt$value,opt$counts[1],opt$convergence)
        }
        
        
        return(results)
    }
    
    
    
    
    

    
    
    
    
    
    
    
    
    
    
    